/*=========================================================================

Program:   Virtual Surgery
Language:  C++
Date:      2009-9-1
Version:   0.1
Mail:		   nelfylee@gmail.com

Copyright (c) 2009 Nelfy Lee. All rights reserved.

=========================================================================*/
#ifndef _SEGMENTATION_THREADS_H_
#define _SEGMENTATION_THREADS_H_

#include "Common.h"
#include "VolumeSource.h"
#include "FilteringThreadBase.h"
#include "itkConnectedThresholdImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkThresholdSegmentationLevelSetImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkImageFileWriter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkShrinkImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkTimeProbe.h"
#include "itkResampleImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkImageRegionIterator.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkConfidenceConnectedImageFilter.h"
#include "itkInvertIntensityImageFilter.h"
#include "itkBinaryContourImageFilter.h"
#include "itkBinaryPruningImageFilter.h"
#include "itkOtsuThresholdImageFilter.h"
#include "itkHessian3DToVesselnessMeasureImageFilter.h"
#include "itkMultiScaleHessianBasedMeasureImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"
#include "itkImageMaskSpatialObject.h"
#include "itkHessianToObjectnessMeasureImageFilter.h"

#include <QList>
#include <QTest>


template<class TImageType>
class Slice
{
public:
	typename TImageType::Pointer origImage;
	ImageUChar2D::Pointer segmentedImage;
	typename TImageType::Pointer contourImage;

	typename TImageType::IndexType seed;
};


enum ThresholdingMethods
{
	thNONE = 0,
	thConnectedThresholding,						// connected threshold
	thConfidenceConnected,						// confidence connected	
	thOtsuThresholding,
	thHessianProcessing,
	thExtractMAXROI
};


template<class CInputPixelType>
class RegionGrowingThread : public FilteringThreadBase
{	
	typedef itk::Image<CInputPixelType, 2> CInputImageType2D;
	typedef itk::Image<CInputPixelType, 3> CInputImageType;


	typedef itk::ConnectedThresholdImageFilter<CInputImageType, CInputImageType> CTFilterType;
	typedef itk::ConfidenceConnectedImageFilter<CInputImageType, CInputImageType> CCFilterType;
	typedef itk::OtsuThresholdImageFilter<CInputImageType, CInputImageType> OtsuFilterType;

public:
	RegionGrowingThread(Source *source, ThresholdingMethods method = thNONE):FilteringThreadBase(source)
	{
		m_lower = 0;
		m_upper = 0;
		m_ccMultiplier = 0;
		m_ccIteration = 0;
		m_method = method;
		m_replaceValue = 1;

		m_inputImage = ((VolumeSource<CInputImageType>*)m_source)->getVolume();
	}
	~RegionGrowingThread()
	{
	}

	virtual void addSeed(typename CInputImageType::IndexType seed)
	{
		m_seedList.push_back(seed);
	}

	virtual void setThreshold(int lower, int upper)
	{
		m_lower = lower;
		m_upper = upper;
	}

	virtual void setReplaceValue(int value = 1)
	{
		m_replaceValue = value;
	}

	virtual void setCCPara(double multiplier, int iteration)
	{
		m_ccMultiplier = multiplier;
		m_ccIteration = iteration;
	}

	void SetHessianParameters(double sigma, double a1, double a2)
	{
		this->sigma = sigma;
		this->a1 = a1;
		this->a2 = a2;
	}

	double sigma;
	double a1;
	double a2;

protected:

	virtual void run()
	{
		typename CInputImageType::Pointer segmented;
		qDebug("[RegionGrowingThread::run]: method is %d", m_method);



		if ( m_method == thConnectedThresholding)
		{
			if ( m_seedList.isEmpty()) 
			{
				qDebug("[Warning]: Select a seed first.");
				return;
			}

			if ( m_lower > m_upper)
			{
				qDebug("Threshold is error.");
				return;
			}

			segmented = connectedThresholdFiltering();
		}
		else if ( m_method == thConfidenceConnected)
		{
			if ( m_seedList.isEmpty()) 
			{
				qDebug("[Warning]: Select a seed first.");
				return;
			}

			if ( m_ccMultiplier == 0 || m_ccIteration == 0)
			{
				qDebug("Confidence Connected parameters didn't set");
				return;
			}

			segmented = confidenceConnectedFiltering();
		}
		else if ( m_method == thOtsuThresholding)
		{
			segmented = otsuFiltering();
		}
		else if ( m_method == thHessianProcessing)
		{
			segmented = hessianFiltering();
		}
		else if ( m_method == thExtractMAXROI)
		{
			segmented = extractMaxROI();
		}
		else
		{			
			return;
		}

		if ( m_replaceValue == 0)
		{
			typedef itk::ImageRegionIterator<CInputImageType> ImageRegionIterator;

			ImageRegionIterator segmentedIt(segmented, segmented->GetLargestPossibleRegion());
			ImageRegionIterator sourceIt(m_inputImage, m_inputImage->GetLargestPossibleRegion());

			sourceIt.GoToBegin();segmentedIt.GoToBegin();
			for ( ; !segmentedIt.IsAtEnd(); ++segmentedIt, ++sourceIt)
			{
				if ( segmentedIt.Get())
				{
					segmentedIt.Set(sourceIt.Get());
				}
			}
		}


		VolumeSource<CInputImageType> *subsource = 
			new VolumeSource<CInputImageType>(QString("Segmented#%1").arg(m_source->rowCount()+1));
		subsource->setVolume(segmented);

		m_source->addChild(subsource);
	}

	// combine with confidenceConnectedFiltering
	typename CInputImageType::Pointer connectedThresholdFiltering()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &RegionGrowingThread::ProcessEvent );

		emit emitProgressStart(QString("CT Segmenting..."));

		typename CTFilterType::Pointer filter = CTFilterType::New();
		filter->SetInput(m_inputImage);
		filter->SetSeed(m_seedList.at(0));
		filter->SetLower(m_lower);
		filter->SetUpper(m_upper);

		if ( m_replaceValue != 0)
			filter->SetReplaceValue(m_replaceValue);

		filter->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());
		filter->Update();

		emit emitProgressDone();

		return filter->GetOutput();
	}

	typename CInputImageType::Pointer confidenceConnectedFiltering()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &RegionGrowingThread::ProcessEvent );

		emit emitProgressStart(QString("CC Segmenting..."));

		typename CCFilterType::Pointer filter = CCFilterType::New();
		filter->SetInput(m_inputImage);
		filter->SetSeed(m_seedList.at(0));
		filter->SetMultiplier(m_ccMultiplier);
		filter->SetNumberOfIterations(m_ccIteration);

		if ( m_replaceValue != 0)
			filter->SetReplaceValue(m_replaceValue);

		filter->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());
		filter->Update();

		emit emitProgressDone();
		return filter->GetOutput();
	}

	typename CInputImageType::Pointer otsuFiltering()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &RegionGrowingThread::ProcessEvent );

		emit emitProgressStart(QString("Otsu Segmenting..."));

		typename OtsuFilterType::Pointer filter = OtsuFilterType::New();
		filter->SetInput(m_inputImage);

		filter->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());
		filter->Update();

		emit emitProgressDone();
		return filter->GetOutput();
	}


	typename CInputImageType::Pointer extractMaxROI()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &RegionGrowingThread::ProcessEvent );

		emit emitProgressStart(QString("Extract maximum roi..."));

		typedef itk::CastImageFilter<CInputImageType, ImageUChar3D> CastImageFilter_Input2UChar3D;
		
		typedef itk::RegionOfInterestImageFilter<CInputImageType, CInputImageType> RegionOfInterestImageFilter;
		typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;

		CastImageFilter_Input2UChar3D::Pointer cast_uchar3d = CastImageFilter_Input2UChar3D::New();
		cast_uchar3d->SetInput(m_inputImage);
		cast_uchar3d->Update();

		typename ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();
		maskSO->SetImage(cast_uchar3d->GetOutput());
		maskSO->Update();

		typename CInputImageType::RegionType regionType = maskSO->GetAxisAlignedBoundingBoxRegion();
 		typename RegionOfInterestImageFilter::Pointer roiFilter = RegionOfInterestImageFilter::New();
		roiFilter->SetInput(m_inputImage);
		roiFilter->SetRegionOfInterest(regionType);
		roiFilter->Update();

		emit emitProgressDone();
		return roiFilter->GetOutput();
	}



	typename CInputImageType::Pointer hessianFiltering2()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &RegionGrowingThread::ProcessEvent );

		emit emitProgressStart(QString("Hessian Segmenting..."));

		typedef itk::Image<float,3> FloatImage;
		typedef itk::ImageFileWriter<FloatImage> FileWriterType;
		typedef itk::CastImageFilter<CInputImageType, FloatImage> CastImageFilter_Input2Float;
		typedef itk::SymmetricSecondRankTensor<float, 3> HessianPixelType;
		typedef itk::Image<HessianPixelType, 3> HessianImage;
		typedef itk::HessianToObjectnessMeasureImageFilter<HessianImage, FloatImage> SatoVesselness;
		typedef itk::MultiScaleHessianBasedMeasureImageFilter<FloatImage, HessianImage, FloatImage> MultiScaleEnhancementFilterType;
		typedef itk::CastImageFilter<FloatImage,CInputImageType> CastImageFilter_Float2Input;

		CastImageFilter_Input2Float::Pointer cast_Input2Float = CastImageFilter_Input2Float::New();
		cast_Input2Float->SetInput(m_inputImage);
		cast_Input2Float->Update();

		SatoVesselness::Pointer satovesselness = SatoVesselness::New();
		satovesselness->SetScaleObjectnessMeasure(false);
		satovesselness->SetBrightObject(true);
		satovesselness->SetAlpha(0.5);
		satovesselness->SetBeta(0.5);	
		satovesselness->SetGamma(10.0);

		MultiScaleEnhancementFilterType::Pointer multiScaleEnhancementFilter = MultiScaleEnhancementFilterType::New();
 		multiScaleEnhancementFilter->SetInput(cast_Input2Float->GetOutput());
		multiScaleEnhancementFilter->SetHessianToMeasureFilter(satovesselness);
		//multiScaleEnhancementFilter->SetSigmaStepMethodToLogarithmic();
		multiScaleEnhancementFilter->SetNumberOfSigmaSteps(10);

 		multiScaleEnhancementFilter->SetSigmaMinimum(0.5);
 		multiScaleEnhancementFilter->SetSigmaMaximum(4.0);
 		//multiScaleEnhancementFilter->SetNumberOfSigmaSteps(10);
 		//multiScaleEnhancementFilter->SetGenerateScalesOutput(true);
		//multiScaleEnhancementFilter->SetGenerateHessianOutput( true );
		multiScaleEnhancementFilter->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

		try
		{
			multiScaleEnhancementFilter->Update();
		}
		catch (itk::ExceptionObject &e)
		{
			std::cerr << e << std::endl;
		}

/*
		FileWriterType::Pointer writer = FileWriterType::New();

		writer->SetFileName("E:/hessian_output.vtk");
		writer->UseCompressionOn();
		writer->SetInput( multiScaleEnhancementFilter->GetOutput() );
		writer->Update();

		writer->SetFileName("E:/hessian_scales_output.vtk");
		writer->UseCompressionOn();
		writer->SetInput( multiScaleEnhancementFilter->GetScalesOutput() );
		writer->Update();

		multiScaleEnhancementFilter->SetSigmaStepMethod( 
			MultiScaleEnhancementFilterType::EquispacedSigmaSteps);
		multiScaleEnhancementFilter->Update();

		FileWriterType::Pointer writer2 = FileWriterType::New();
		writer2->SetFileName("E:/hessian_equispace.vtk");
		writer2->UseCompressionOn();
		writer2->SetInput( multiScaleEnhancementFilter->GetOutput() );
		writer2->Update();

		multiScaleEnhancementFilter->NonNegativeHessianBasedMeasureOff(); 
		multiScaleEnhancementFilter->Update();

		FileWriterType::Pointer writer3 = FileWriterType::New();
		writer3->SetFileName("E:/hessian_nonnegative.vtk");
		writer3->UseCompressionOn();
		writer3->SetInput( multiScaleEnhancementFilter->GetScalesOutput() );
		writer3->Update();
		*/

		
		typedef itk::RescaleIntensityImageFilter<FloatImage, FloatImage> RescaleIntensityImageFilter;
		RescaleIntensityImageFilter::Pointer rescale = RescaleIntensityImageFilter::New();
		rescale->SetInput(multiScaleEnhancementFilter->GetOutput());
		rescale->SetOutputMaximum(255);
		rescale->SetOutputMinimum(0);

 		CastImageFilter_Float2Input::Pointer cast_Float2Input = CastImageFilter_Float2Input::New();
 		cast_Float2Input->SetInput(rescale->GetOutput());
 		cast_Float2Input->Update();

		emit emitProgressDone();
		return cast_Float2Input->GetOutput();
	}




	typename CInputImageType::Pointer hessianFiltering()
	{
		ProgressCommandType::Pointer progressCmd = ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &RegionGrowingThread::ProcessEvent );

		emit emitProgressStart(QString("Hessian Segmenting..."));

		typedef itk::Image<float,3> FloatImage;
		typedef itk::Image<double,3> DoubleImage;
		typedef itk::CastImageFilter<CInputImageType, DoubleImage> CastImageFilter;
		//typedef itk::SymmetricSecondRankTensor<double, 3> HessianPixelType;
		typedef itk::Image<double, 3> HessianImage;
		typedef itk::Hessian3DToVesselnessMeasureImageFilter<double> SatoVesselness;
		typedef itk::MultiScaleHessianBasedMeasureImageFilter<FloatImage, HessianImage, FloatImage> MultiScaleEnhancementFilterType;
		typedef itk::CastImageFilter<DoubleImage,CInputImageType> CastImageFilter2;

		CastImageFilter::Pointer cast1 = CastImageFilter::New();
		cast1->SetInput(m_inputImage);
		cast1->Update();

		typedef   itk::HessianRecursiveGaussianImageFilter<DoubleImage > HessianFilterType;
		HessianFilterType::Pointer hessian = HessianFilterType::New();
		hessian->SetSigma( this->sigma );
		hessian->SetInput( cast1->GetOutput() );
		hessian->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

		typedef   itk::Hessian3DToVesselnessMeasureImageFilter< 
			double > VesselnessMeasureFilterType;
		VesselnessMeasureFilterType::Pointer vesselness = 
			VesselnessMeasureFilterType::New();
		vesselness->SetInput( hessian->GetOutput() );
		vesselness->SetAlpha1(this->a1);
		vesselness->SetAlpha2(this->a2);
		vesselness->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

		CastImageFilter2::Pointer cast2 = CastImageFilter2::New();
		cast2->SetInput(vesselness->GetOutput());
		cast2->Update();

		emit emitProgressDone();
		return cast2->GetOutput();
	}


protected:
	ThresholdingMethods m_method;
	int m_replaceValue;
	double m_ccMultiplier;
	int m_ccIteration;
	
	QList<typename CInputImageType::IndexType> m_seedList;

	typename CInputImageType::Pointer m_inputImage;
	typename CTFilterType::InputImagePixelType m_lower;
	typename CTFilterType::InputImagePixelType m_upper;
};


enum LevelSetMode
{
	lsm2DSlice = 2,
	lsm3DVolume
};


typedef struct _ThresholdLevelsetPara
{
	float tlsDist;
	float tlsRMS;
	float tlsCurv;
	float tlsISO;
	int tlsIteration;
	float tlsProp;
	int tlsGoLeft;
	int tlsGoRight;
}TLS_PARA;



template<class TPixelType>
class ThresholdLevelSetThread : public RegionGrowingThread<TPixelType>
{
	typedef itk::Image<TPixelType, 3> TImageType;
	typedef itk::Image<TPixelType, 2> TImageType2D;

	typedef itk::ShrinkImageFilter<TImageType, TImageType> ShrinkType;
	typedef itk::CastImageFilter<TImageType, ImageFloat3D> CastToRealFilterType;
	typedef  itk::ThresholdSegmentationLevelSetImageFilter< TImageType2D, 
		ImageFloat2D> ThresholdSegmentationLevelSetImageFilterType2;
	typedef  itk::ThresholdSegmentationLevelSetImageFilter< TImageType, 
		ImageFloat3D> ThresholdSegmentationLevelSetImageFilterType3D;

	
	
	typedef itk::BinaryThresholdImageFilter<ImageFloat3D, TImageType>	ThresholdingFilterType3D;
	typedef itk::ImageRegionIterator<TImageType2D> ImageRegionIterator;
	typedef itk::ImageRegionIterator<TImageType> ImageRegionIterator3D;

	typedef itk::CurvatureFlowImageFilter<TImageType2D, ImageFloat2D> CurvatureFlowImageFilter;
	typedef itk::CurvatureFlowImageFilter<TImageType, ImageFloat3D> CurvatureFlowImageFilter3D;

public:
	ThresholdLevelSetThread(Source *source):RegionGrowingThread<TPixelType>(source){}
	~ThresholdLevelSetThread(){}

	// Choose 3D volume segmentation or 2D slicing segmentation
	virtual void setMode(LevelSetMode mode){ m_mode = mode;}

	// Set Level set parameters
	virtual void setParameters(TLS_PARA& para) { memcpy(&m_para, &para, sizeof(TLS_PARA));}

protected:
	virtual void run()
	{
		m_sourceVolume = ((VolumeSource<TImageType>*)RegionGrowingThread<TPixelType>::m_source)->getVolume();

		if (m_mode == lsm2DSlice)
		{
			emit FilteringThreadBase::emitProgressStart("Liver 2D slicing segmentation...");
			run2D();
		}
		else if ( m_mode == lsm3DVolume)
		{
			emit FilteringThreadBase::emitProgressStart("Liver 3D slicing segmentation...");
			run3D();
			emit FilteringThreadBase::emitProgressDone();
		}
		else
		{
			qDebug("[ThresholdLevelSetThread]: Mode error");
			return;
		}

		emit FilteringThreadBase::emitProgressDone();
	}


private:
	void run2D()
	{
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typedef itk::CastImageFilter<TImageType, ImageUChar3D> CastImageFilter;
		typedef itk::ImageRegionIterator<ImageUChar3D> ImageRegionIterator;
		typedef itk::ImageRegionIterator<ImageUChar2D> ImageRegionIterator2D;
		typedef itk::ImageToVTKImageFilter<ImageUChar2D> ConnectorUChar2D;
		
		typename ImageDuplicator::Pointer duplicator = ImageDuplicator::New();
		typename CastImageFilter::Pointer caster = CastImageFilter::New();
		
		duplicator->SetInputImage(m_sourceVolume);
		duplicator->Update();
		duplicator->GetOutput()->FillBuffer(0);

		caster->SetInput(duplicator->GetOutput());
		caster->Update();		

		typename ImageUChar3D::Pointer targetVolume = caster->GetOutput();
		typename ImageUChar3D::SpacingType spacing = caster->GetOutput()->GetSpacing();


		VolumeSource<ImageUChar3D> *targetSource = 
        new VolumeSource<ImageUChar3D>(QString("Temp#%1").arg(RegionGrowingThread<TPixelType>::m_source->rowCount()+1));
		targetSource->setVolume(targetVolume);

		RegionGrowingThread<TPixelType>::m_source->addChild(targetSource);


		emit FilteringThreadBase::emitStartRealTimeSlicing(targetSource);

		typename TImageType::RegionType inputRegion = m_sourceVolume->GetLargestPossibleRegion();
		typename TImageType::SizeType  inputSize = inputRegion.GetSize();

		//int totalSlices = inputSize[2];
		int sliceNum = RegionGrowingThread<TPixelType>::m_seedList.at(0)[2];


		segmentSlice(m_sourceVolume, targetVolume, sliceNum, 0);
		ImageUChar2D::Pointer targetSlice = extractSlice<ImageUChar3D, ImageUChar2D>(targetVolume, sliceNum);

		ConnectorUChar2D::Pointer connector = ConnectorUChar2D::New();
		connector->SetInput(targetSlice);
		connector->Update();

		vtkImageData *m_vtkImage = vtkImageData::New();
		m_vtkImage->DeepCopy(connector->GetOutput());

		emit FilteringThreadBase::emitUpdateSlicing(m_vtkImage, sliceNum);

		bool validSlice = false;

		if ( m_para.tlsGoRight != -1)
		{
			//while( sliceNum++ < totalSlices - 1)
			while( sliceNum++ < m_para.tlsGoRight)
			{
				validSlice = false;

				segmentSlice(m_sourceVolume, targetVolume, sliceNum, 1);

				targetSlice = extractSlice<ImageUChar3D, ImageUChar2D>(targetVolume, sliceNum);
				ImageRegionIterator2D segmentedIt(targetSlice, targetSlice->GetLargestPossibleRegion());
				for ( segmentedIt.GoToBegin(); !segmentedIt.IsAtEnd(); ++segmentedIt)
				{
					if ( segmentedIt.Get())
					{
						validSlice = true;
						break;
					}
				}

				if ( validSlice == false) 
				{
					qDebug() << "Skip slice " << sliceNum;
					break;
				}

				ConnectorUChar2D::Pointer connector = ConnectorUChar2D::New();
				connector->SetInput(targetSlice);
				connector->Update();

				vtkImageData *m_vtkImage = vtkImageData::New();
				m_vtkImage->DeepCopy(connector->GetOutput());

				emit FilteringThreadBase::emitUpdateSlicing(m_vtkImage, sliceNum);
			}
		}

		sliceNum = RegionGrowingThread<TPixelType>::m_seedList.at(0)[2];

		if ( m_para.tlsGoLeft != -1)
		{
			//while( sliceNum-- > 0)
			while( sliceNum-- > m_para.tlsGoLeft)
			{
				validSlice=false;

				segmentSlice(m_sourceVolume, targetVolume, sliceNum, -1);

				targetSlice = extractSlice<ImageUChar3D, ImageUChar2D>(targetVolume, sliceNum);
				ImageRegionIterator2D segmentedIt(targetSlice, targetSlice->GetLargestPossibleRegion());
				for ( segmentedIt.GoToBegin(); !segmentedIt.IsAtEnd(); ++segmentedIt)
				{
					if ( segmentedIt.Get())
					{
						validSlice = true;
						break;
					}
				}

				if ( validSlice == false) 
				{
					qDebug() << "Skip slice " << sliceNum;
					break;
				}

				ConnectorUChar2D::Pointer connector = ConnectorUChar2D::New();
				connector->SetInput(targetSlice);
				connector->Update();

				vtkImageData *m_vtkImage = vtkImageData::New();
				m_vtkImage->DeepCopy(connector->GetOutput());

				emit FilteringThreadBase::emitUpdateSlicing(m_vtkImage, sliceNum);
			}
		}

		targetSource->setVolume(targetVolume);

		emit FilteringThreadBase::emitStopRealTimeSlicing();
		return;
	}


	template<class FInputImageType, class FOutputImageType>
	typename FOutputImageType::Pointer extractSlice(typename FInputImageType::Pointer volume, int num)
	{
		typedef itk::ExtractImageFilter<FInputImageType, FOutputImageType> SliceExtractor;
		typename SliceExtractor::Pointer sliceExtractor = SliceExtractor::New();


		typename FInputImageType::RegionType inputRegion = volume->GetLargestPossibleRegion();
		typename FInputImageType::SizeType  inputSize = inputRegion.GetSize();
		inputSize[2] = 0;


		typename FInputImageType::IndexType start = inputRegion.GetIndex();
		start[2] = num;

		typename FInputImageType::RegionType desiredRegion;
		desiredRegion.SetSize(inputSize);
		desiredRegion.SetIndex(start);

		sliceExtractor->SetExtractionRegion(desiredRegion);
		sliceExtractor->SetInput(volume);
		sliceExtractor->Update();

		return sliceExtractor->GetOutput();
	}


	template<class TInputImageType, class TOutputImageType>
	typename TOutputImageType::PixelType* getSlicePointer(typename TInputImageType::Pointer volume, int num)
	{
		typename TInputImageType::RegionType inputRegion = volume->GetLargestPossibleRegion();
		typename TInputImageType::SizeType  inputSize = inputRegion.GetSize();
		inputSize[2] = 0;


		typename TInputImageType::PixelType *pointer = volume->GetBufferPointer();

		return (pointer + inputSize[0]*inputSize[1]*num);
	}


	void fillHoles(typename ImageUChar2D::Pointer reference, typename ImageUChar2D::PixelType *target)
	{
		typedef itk::InvertIntensityImageFilter<ImageUChar2D> InvertIntensityImageFilter;
		typedef itk::ConnectedThresholdImageFilter<ImageUChar2D, ImageUChar2D> ConnectedThresholdImageFilter;

		InvertIntensityImageFilter::Pointer invert = InvertIntensityImageFilter::New();
		ConnectedThresholdImageFilter::Pointer filter = ConnectedThresholdImageFilter::New();

		invert->SetInput(reference);
		invert->SetMaximum(1);
		
		filter->SetInput(invert->GetOutput());

		static ImageUChar2D::IndexType seed2;
		seed2[0] = 1;
		seed2[1] = 1;

		filter->SetSeed(seed2);
		filter->SetLower(1);
		filter->SetUpper(1);
		filter->Update();

		invert->SetInput(filter->GetOutput());
		invert->Update();

		ImageUChar2D::RegionType outputRegion = invert->GetOutput()->GetLargestPossibleRegion();
		ImageUChar2D::SizeType  outputSize = outputRegion.GetSize();

		qMemCopy(target, invert->GetOutput()->GetBufferPointer(),outputSize[0]*outputSize[1]);
	}


	void segmentSlice(typename TImageType::Pointer sourceVolume, 
		typename ImageUChar3D::Pointer targetVolume, 
		int sliceNum, 
		int direction						// -1, 0, 1
		)
	{
		typename TImageType2D::Pointer sourceSlice = extractSlice<TImageType,TImageType2D>(sourceVolume, sliceNum);
		typename ImageUChar2D::PixelType *targetSlice = getSlicePointer<ImageUChar3D,ImageUChar2D>(targetVolume, sliceNum);
		typename TImageType2D::Pointer maskSlice = extractSlice<ImageUChar3D,TImageType2D>(targetVolume, sliceNum-direction);

		if ( direction == 0)
		{
			typename TImageType2D::IndexType seed;
			
			seed[0] = RegionGrowingThread<TPixelType>::m_seedList.at(0)[0];
			seed[1] = RegionGrowingThread<TPixelType>::m_seedList.at(0)[1];

			segmentSliceBySeed<TImageType2D, ImageUChar2D, 2>(sourceSlice,targetSlice,	seed);

			return;
		}
		else
		{
			segmentSliceByMask(sourceSlice, targetSlice, maskSlice);
		}


		ImageUChar2D::Pointer reference = extractSlice<ImageUChar3D,ImageUChar2D>(targetVolume, sliceNum);

		fillHoles(reference, targetSlice);

		qDebug("Segment slice %d, reference slice %d", sliceNum, sliceNum-direction);
		qDebug("Done");
	}


	void run3D()
	{
		typename TImageType::IndexType seed;

		seed[0] = RegionGrowingThread<TPixelType>::m_seedList.at(0)[0];
		seed[1] = RegionGrowingThread<TPixelType>::m_seedList.at(0)[1];
		seed[2] = RegionGrowingThread<TPixelType>::m_seedList.at(0)[2];

		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typedef itk::CastImageFilter<TImageType, ImageUChar3D> CastImageFilter;

		typename ImageDuplicator::Pointer duplicator = ImageDuplicator::New();
		typename CastImageFilter::Pointer caster = CastImageFilter::New();

		duplicator->SetInputImage(m_sourceVolume);
		duplicator->Update();
		duplicator->GetOutput()->FillBuffer(0);

		caster->SetInput(duplicator->GetOutput());
		caster->Update();

		ImageUChar3D::PixelType *targetVolume = caster->GetOutput()->GetBufferPointer();

		segmentSliceBySeed<TImageType, ImageUChar3D, 3>(m_sourceVolume,targetVolume,	seed);

		VolumeSource<ImageUChar3D> *targetSource = 
			new VolumeSource<ImageUChar3D>(QString("Temp#%1").arg(RegionGrowingThread<TPixelType>::m_source->rowCount()+1));
		targetSource->setVolume(caster->GetOutput());


		RegionGrowingThread<TPixelType>::m_source->addChild(targetSource);

	}





	int segmentSliceByMask(typename TImageType2D::Pointer source,
		typename ImageUChar2D::PixelType *target,
		typename TImageType2D::Pointer mask)
	{
		typedef itk::CastImageFilter<TImageType2D, ImageFloat2D> CastImageFilter;
		typedef itk::BinaryThresholdImageFilter<ImageFloat2D, ImageUChar2D>	ThresholdingFilterType;

		static double timeTaken;
		static itk::TimeProbe smoothClock, segClock;

		segClock.Start();	qDebug("Segmenting...");

		typename RegionGrowingThread<TPixelType>::ProgressCommandType::Pointer progressCmd = RegionGrowingThread<TPixelType>::ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &ThresholdLevelSetThread::ProcessEvent );


		ImageRegionIterator maskIt(mask, mask->GetRequestedRegion());
		
		for ( maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
		{
			if ( maskIt.Get())
				maskIt.Set(0);
			else
				maskIt.Set(1);
		}

		typename CastImageFilter::Pointer caster = CastImageFilter::New();
		caster->SetInput(source);

		ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();

		thresholder->SetLowerThreshold( -1000.0 );
		thresholder->SetUpperThreshold(     0.0 );

		thresholder->SetOutsideValue(  0  );
		thresholder->SetInsideValue(  1 );

		typename ThresholdSegmentationLevelSetImageFilterType2::Pointer thresholdSegmentation =
			ThresholdSegmentationLevelSetImageFilterType2::New();

		thresholdSegmentation->SetPropagationScaling( m_para.tlsProp );
		thresholdSegmentation->SetCurvatureScaling(m_para.tlsCurv);

		thresholdSegmentation->SetMaximumRMSError( m_para.tlsRMS );
		thresholdSegmentation->SetNumberOfIterations( 800 );

		thresholdSegmentation->SetLowerThreshold( (double)RegionGrowingThread<TPixelType>::m_lower );
		thresholdSegmentation->SetUpperThreshold( (double)RegionGrowingThread<TPixelType>::m_upper);
		thresholdSegmentation->SetIsoSurfaceValue(0.5);


		thresholdSegmentation->SetInput( mask );
		thresholdSegmentation->SetFeatureImage( caster->GetOutput());
		thresholder->SetInput( thresholdSegmentation->GetOutput() );

		thresholdSegmentation->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

		try
		{
			thresholder->Update();
		}
		catch ( itk::ExceptionObject & excep )
		{
			std::cerr << "Exception caught !" << std::endl;
			std::cerr << excep << std::endl;
		}

		std::cout << std::endl;
		std::cout << "No. elpased iterations: " << thresholdSegmentation->GetElapsedIterations() << std::endl;
		std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() << std::endl;


		segClock.Stop();
		timeTaken = segClock.GetMeanTime();
		std::cout << "[OK]:" << std::endl;
		std::cout << "[ Time ] Segmentation took " << timeTaken << " s." << std::endl;

		typename TImageType2D::RegionType inputRegion = thresholder->GetOutput()->GetLargestPossibleRegion();
		typename TImageType2D::SizeType  inputSize = inputRegion.GetSize();

		qMemCopy(target, thresholder->GetOutput()->GetBufferPointer(), inputSize[0]*inputSize[1]);

		return 0;
	}



	int segmentSliceByMask(typename TImageType::Pointer source,
		//typename TImageType::PixelType *target,
		typename TImageType::Pointer mask)
	{
		static double timeTaken;
		static itk::TimeProbe smoothClock, segClock;

		segClock.Start();	qDebug("Final matching...");

		typename RegionGrowingThread<TPixelType>::ProgressCommandType::Pointer progressCmd = RegionGrowingThread<TPixelType>::ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &ThresholdLevelSetThread::ProcessEvent );


		ImageRegionIterator3D maskIt(mask, mask->GetLargestPossibleRegion());

		for ( maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt)
		{
			if ( maskIt.Get())
				maskIt.Set(0);
			else
				maskIt.Set(1);
		}

		typename CurvatureFlowImageFilter3D::Pointer smoothing = CurvatureFlowImageFilter3D::New();
		smoothing->SetInput(source);
		smoothing->SetNumberOfIterations(0);
		smoothing->SetTimeStep(0.0625);



		typename ThresholdingFilterType3D::Pointer thresholder = ThresholdingFilterType3D::New();

		thresholder->SetLowerThreshold( -1000.0 );
		thresholder->SetUpperThreshold(     0.0 );

		thresholder->SetOutsideValue(  0  );
		thresholder->SetInsideValue(  1 );

		typename ThresholdSegmentationLevelSetImageFilterType3D::Pointer thresholdSegmentation =
			ThresholdSegmentationLevelSetImageFilterType3D::New();

		thresholdSegmentation->SetPropagationScaling( m_para.tlsProp );
		thresholdSegmentation->SetCurvatureScaling(m_para.tlsCurv);

		thresholdSegmentation->SetMaximumRMSError( m_para.tlsRMS );
		thresholdSegmentation->SetNumberOfIterations( 1200 );

		thresholdSegmentation->SetLowerThreshold( (double)RegionGrowingThread<TPixelType>::m_lower );
		thresholdSegmentation->SetUpperThreshold( (double)RegionGrowingThread<TPixelType>::m_upper);
		thresholdSegmentation->SetIsoSurfaceValue(0.5);


		thresholdSegmentation->SetInput( mask );
		thresholdSegmentation->SetFeatureImage( smoothing->GetOutput());
		thresholder->SetInput( thresholdSegmentation->GetOutput() );

		thresholdSegmentation->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

		try
		{
			thresholder->Update();
		}
		catch ( itk::ExceptionObject & excep )
		{
			std::cerr << "Exception caught !" << std::endl;
			std::cerr << excep << std::endl;
		}

		std::cout << std::endl;
		std::cout << "No. elpased iterations: " << thresholdSegmentation->GetElapsedIterations() << std::endl;
		std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() << std::endl;


		segClock.Stop();
		timeTaken = segClock.GetMeanTime();
		std::cout << "[OK]:" << std::endl;
		std::cout << "[ Time ] Segmentation took " << timeTaken << " s." << std::endl;

		typename TImageType::RegionType inputRegion = thresholder->GetOutput()->GetLargestPossibleRegion();
		typename TImageType::SizeType  inputSize = inputRegion.GetSize();

		typedef itk::ImageFileWriter<TImageType> ImageFileWriter;
		typename ImageFileWriter::Pointer writer = ImageFileWriter::New();
		writer->SetInput(thresholder->GetOutput());
		writer->SetFileName("matched.vtk");

		writer->Update();

//		qMemCopy(target, thresholder->GetOutput()->GetBufferPointer(), inputSize[0]*inputSize[1]*sizeof(TPixelType));

		return 0;
	}




	template<class FInputImage, class FOutputImage, unsigned int Dim>
	int segmentSliceBySeed(typename FInputImage::Pointer inputImage,
		typename FOutputImage::PixelType *outputImage,
		typename FInputImage::IndexType seed)
	{
 		typedef itk::Image<float, Dim> InternalImageType;

		typedef itk::CastImageFilter<FInputImage, InternalImageType> CastImageFilter;
		typedef itk::ThresholdSegmentationLevelSetImageFilter< InternalImageType, 
			InternalImageType> ThresholdSegmentationLevelSetImageFilterType;
		typedef itk::BinaryThresholdImageFilter<InternalImageType, FOutputImage>	ThresholdingFilterType;
		typedef  itk::FastMarchingImageFilter<InternalImageType,InternalImageType> FastMarchingFilterType;


		double timeTaken;
		itk::TimeProbe smoothClock, segClock;

		segClock.Start();

		typename RegionGrowingThread<TPixelType>::ProgressCommandType::Pointer progressCmd = RegionGrowingThread<TPixelType>::ProgressCommandType::New();
		progressCmd->SetCallbackFunction( this, &ThresholdLevelSetThread::ProcessEvent );

		typename CastImageFilter::Pointer caster = CastImageFilter::New();
		caster->SetInput(inputImage);

		typename ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();

		thresholder->SetLowerThreshold( -1000.0 );
		thresholder->SetUpperThreshold(     0.0 );

		thresholder->SetOutsideValue(0);
		thresholder->SetInsideValue(1);

		typename FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
		typename ThresholdSegmentationLevelSetImageFilterType::Pointer thresholdSegmentation =
			ThresholdSegmentationLevelSetImageFilterType::New();

		thresholdSegmentation->SetPropagationScaling( m_para.tlsProp );
		thresholdSegmentation->SetCurvatureScaling(m_para.tlsCurv);

		thresholdSegmentation->SetMaximumRMSError( m_para.tlsRMS );
		thresholdSegmentation->SetNumberOfIterations( m_para.tlsIteration );

		thresholdSegmentation->SetLowerThreshold( (double)RegionGrowingThread<TPixelType>::m_lower );
		thresholdSegmentation->SetUpperThreshold( (double)RegionGrowingThread<TPixelType>::m_upper);
		thresholdSegmentation->SetIsoSurfaceValue(0.5);


		thresholdSegmentation->SetInput( fastMarching->GetOutput() );
		thresholdSegmentation->SetFeatureImage( caster->GetOutput());
		thresholder->SetInput( thresholdSegmentation->GetOutput() );



		typedef typename FastMarchingFilterType::NodeContainer           NodeContainer;
		typedef typename FastMarchingFilterType::NodeType                NodeType;
        typename NodeContainer::Pointer seeds = NodeContainer::New();

		typename InternalImageType::IndexType  seedPosition;

		seedPosition[0] = seed[0];
		seedPosition[1] = seed[1];
		if ( Dim == 3)
			seedPosition[2] = seed[2];

		NodeType node;
		const double seedValue = - m_para.tlsDist;

		node.SetValue( seedValue );
		node.SetIndex( seedPosition );

		seeds->Initialize();
		seeds->InsertElement( 0, node );

		fastMarching->SetTrialPoints(  seeds  );
		fastMarching->SetSpeedConstant( 1.0 );

		fastMarching->SetOutputSize(inputImage->GetBufferedRegion().GetSize() );
		fastMarching->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

// 		if ( Dim == 2)
// 			fastMarching->SetStoppingValue(100);

		thresholdSegmentation->AddObserver( itk::ProgressEvent(), progressCmd.GetPointer());

		try
		{
			thresholder->Update();
		}
		catch ( itk::ExceptionObject & excep )
		{
			std::cerr << "Exception caught !" << std::endl;
			std::cerr << excep << std::endl;
		}

/*
		typedef itk::ImageFileWriter<InternalImageType> ImageFileWriter;
		ImageFileWriter::Pointer writer = ImageFileWriter::New();

		writer->SetInput(fastMarching->GetOutput());
		writer->SetFileName("../fastMarching.vtk");
		writer->Update();

		ImageFileWriter::Pointer writer2 = ImageFileWriter::New();
		writer2->SetInput(thresholdSegmentation->GetOutput());
		writer2->SetFileName("../thresholdSegmentation.vtk");
		writer2->Update();

		typedef itk::ImageFileWriter<FOutputImage> ImageFileWriter2;
		ImageFileWriter2::Pointer writer3 = ImageFileWriter2::New();
		writer3->SetInput(thresholder->GetOutput());
		writer3->SetFileName("../thresholder.vtk");
		writer3->Update();
*/


		qDebug("No. elpased iterations: %d", thresholdSegmentation->GetElapsedIterations());
		qDebug("RMS change: %f", thresholdSegmentation->GetRMSChange());

		segClock.Stop();
		timeTaken = segClock.GetMeanTime();
		std::cout << "Segmentation took " << timeTaken << " s." << std::endl;

		typename FOutputImage::RegionType outputRegion = thresholder->GetOutput()->GetLargestPossibleRegion();
		typename FOutputImage::SizeType  outputSize = outputRegion.GetSize();

		if ( Dim == 2)
		{
			qMemCopy(outputImage, thresholder->GetOutput()->GetBufferPointer(), 
				outputSize[0]*outputSize[1]*sizeof(typename FOutputImage::PixelType));
		}
		else
		{
			qMemCopy(outputImage, thresholder->GetOutput()->GetBufferPointer(), 
				outputSize[0]*outputSize[1]*outputSize[2]*sizeof(typename FOutputImage::PixelType));
		}

		return 0;
}








private:
	LevelSetMode m_mode;
	TLS_PARA m_para;
	typename TImageType::Pointer m_sourceVolume;
};






#endif





